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Analytic  Ray  Tracing  Using  Ermakov  Invariants 


1.  INTRODUCTION 

The  impetus  for  developing  a  computational  technique  capable  of  rapidly  producing  ray 
traces  and  lonograma  came  from  the  requirements  of  an  experiment  designed  to  measure  HF 
propagation  in  the  ionospheric  E-F2  duct  in  situ.  In  the  experiment,  a  receiver  placed  on  a 
polar  orbiting  satellite  Is  to  detect  6-30  MHz  signals  injected  into  the  duct  either  by  means  of 
naturally  occurring  ionospheric  gradients  or  by  scattering  from  an  artificially  heated 
ionospheric  volume.  Since  a  time  delay  vs  frequency  ionogram  will  change  depending  on 
position  In  the  orbit,  a  capability  for  producing  essentially  real-time  ionograms  is  needed  to 
direct  the  experiment  efficiently. 

Programs  exist  that  can  numerically  Integrate  the  differential  equations  of  motion  (such 
as  the  Haselgrove  Hamiltonian  equations)  that  describe  a  propagating  ray. 1  Depending  on  the 
path  length,  the  time  required  to  produce  an  ionogram  on  a  mainframe  computer  can  run  from 
20-200  minutes.  The  time  required  to  evaluate  an  ionogram  using  the  analysis  to  follow  will 
require  less  than  100  seconds  and  can  be  run  on  small  computers. 

(Received  for  Publication  14  Nov.  1989) 

1  Jones.  R.M.,  Stephenson,  J.J.  (1975)  A  Versatile  Three-Dimensional  Ray  Tracing  Computer 
Program  for  Radio  Waves  In  the  Ionosphere.  U.S.  OT  Report.  75-76. 
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The  method  developed  here  starts  from  a  description  of  propagating  HF  waves  by  the  scalar 
wave  equation  in  the  paraxial  approximation,  which  is  valid  if  changes  in  the  ionosphere  are 
small  over  distances  on  the  order  of  a  wavelength  (X.  -  10- 50m).  By  redefining  z,  the  distance 
along  the  great  circle  path  as  a  time  coordinate,  the  paraxial  equation  can  be  cast  In  the  form 
of  the  time  dependent  Schroedinger  equation. 

The  time  dependent  Schroedinger  equation  has  recently  been  shown  to  be  soluble  for  a 
certain  class  of  potentials  by  Ermakov  invariant  theory.2  This  class  consists  of  potentials 
obtained  by  canonical  transformations  that  are  parametrized  by  functions  that  obey  a  set  of 
time-dependent  differential  constraints.  By  approximating  a  path-dependent  ionospheric 
potential,  such  that  it  satisfies  these  constraints,  exact  time-dependent  scalar  wave  solutions 
to  the  approximating  ionospheric  potential  may  be  obtained. 

Since  It  is  difficult  to  find  soluble  double-well  potentials  that  can  approximate  the 
ionosphere  over  extended  ranges,  as  a  practical  matter  it  is  necessary  to  Join  two  (or  more) 
single-well  potentials  at  a  common  boundary  to  model  the  total  ionospheric  potential.  A 
consistent  calculation  can  then  be  made  in  the  ray  approximation,  which  corresponds  to 
solving  Newton's  equations  of  motion  for  a  particle  in  a  time-dependent  one-dimensional 
potential  well.  By  approximating  the  ionospheric  potential  along  a  great  circle  path  by  joined 
single-well  potentials,  such  that  each  satisfies  the  differential  constraints  of  Ermakov  theory, 
the  solution  to  the  ray  tracing  problem  is  mapped  onto  the  problem  of  solving  for  the 
equations  of  motion  of  a  particle  in  a  static  potential. 

By  choosing  the  fitting  potential  to  be  a  combination  of  Morse  and  linear  potentials  or 
quadratic  and  linear  potentials  or  any  potential  with  known  solutions  that  when  shifted  and 
scaled  conforms  to  the  ionospheric  potential  over  the  chosen  range,  a  good  approximation  for 
the  entire  ray  tracing  problem  may  thus  be  obtained  analytically.  The  approximations  that 
Inhere  in  this  method  of  tracing  rays  in  a  model  ionosphere  enter  through  the  choice  of  the 
fitting  potential,  the  paraxial  approximation,  and  the  fact  that  the  differential  constraint 
equations  are  only  approximately  obeyed. 

The  outline  of  this  report  is  as  follows:  In  Section  2  we  indicate  how  the  paraxial 
approximation  to  the  scalar  wave  equation  may  be  viewed  as  a  time-dependent  Schroedinger 
equation.  The  use  of  Ermakov  invariant  theory  for  solving  the  time-dependent  Schroedinger 
equation  Is  then  outlined.  In  Section  3  we  discuss  how  the  application  of  the  time-dependent 
Schroedinger  formalism  may  be  applied  to  the  Ionospheric  ray  tracing  problem  by  using 
Morse  and  linear  potentials  to  fit  the  ionosphere.  The  possibility  of  using  double- well 
potentials  obtained  through  the  Darboux  transformation  is  also  discussed.  In  Section  4  we 
discuss  the  relation  of  coherent  states  of  the  time-dependent  Schroedinger  equation  and  rays 
for  the  Morse  potential.  The  possibility  of  multiple  focusing  of  rays  is  pointed  out.  This 
behavior  is  analogous  to  convergence  zones  in  the  case  of  ocean  acoustic  ducts  and  provides  an 
explanation  for  both  cases.  In  Section  5  we  discuss  ray  tracing  in  the  potential  wells  provided 


2  Ray,  J.R.  (1982)  Exact  Solutions  to  the  Time- Dependent  Schroedinger  Equation,  Phys. 
Rev.,  A26:729/Hartley.  J.G.,  Ray,  J.R  (1981)  Ermakov  Systems  and  Quantum-Mechanical 
Superposition  Laws.  Phys.  Reu.,A24:2873. 
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directly  by  IONCAP.3  Examples  of  ray  traces  and  lonograms  calculated  both  by  the  Jones- 
Stephenson  ray  tracing  program  using  IONCAP  and  by  the  present  method  are  given. 


2.  RAT  TRACING  AND  THE  TIME-DEPENDENT  8CHROEDINGER  EQUATION 


Our  description  of  Ionospheric  HF  propagation  starts  from  the  scalar  wave  equation  for 
a  transverse  electric  field  component.  For  propagation  along  a  great  circle  path  this  can  be 
written  as 


d2E 

dr* 


+  i^  + 

r  dr 


xdji 

iW 


n2(r.$)  E(r,$)  *  0 


(1) 


with  k  «=  fil  . 
o  c 


Assuming  a  solution  of  the  form 


E(r.<>) 


E(r,*) 


(2) 


This  can  be  written  as 


d2g  +  L  +  1_ 

dr2  r  <^r  r2 


^  +  2.0^ 

d*2 


+  k2na(r.*)  E(r,*) 
r  *  ° 


0 


(3) 


In  terms  of  the  coordinates  for  height  x  above  the  Earth's  surface  x  =  r  -  Ro  where  Ro  Is  the 
Earth  s  radius,  and  range  z  along  the  great  circle  path  at  the  Earth  s  surface  z  *  Ro$.  this  can  be 
rewritten  as 


3  Lloyd,  J.L.,  Haydon.  G.W.,  Lucas,  D.L..  and  Teters,  L.R.  (1984)  Estimating  the  Performance 
of  Telecommunications  Systems  Using  the  Ionospheric  Transmission  Channel.  Vol.  1. 
Techniques  for  Analyzing  Ionospheric  Effects  Upon  HF  Systems.  NTTA/ITS. 
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(4) 


+ _ L__^  + 

dx2  (R0  +  x)ax 
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2az2  0+x/R0) 


dE 
2  dz 


2  2 

Pk0E  ,* 

- Q - +  k  n  £ 

.  2  0 
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where  the  angular  momentum  P  has  been  re-expressed  as 


P  =  PkoRo‘ 


(5) 


Since  derivatives  In  height  occur  on  a  scale  of  tens  of  km,  the  term  — ^ — -  —  may  be 

(Ro  +  x)  dx 

^2E 

neglected  In  comparison  with  the  term  — —  .  If  we  assume  that  changes  In  n  (x.z)  ((and  E  (x.z)) 

dx~ 

occur  on  a  scale  much  larger  than  the  wavelength  X  (*  2 n/ko)  of  the  propagating  wave,  then  we 
may  neglect  terms  operated  on  by  d/dz  In  comparison  to  those  multiplied  by  lpko.  Expanding 
the  remaining  terms  to  lowest  order  In  x/Ro  yields  the  approximate  wave  equation 


k0  dz  2  k*  dx2 


2  2 
n  (x.  z)  P  x 

2  “  r7 


£ 

2 


e  (x,  z) 


(6) 


which  is  In  the  form  of  the  time-dependent  Schroedlnger  equation  with  k  j1  replacing  h.  z/p 
replacing  the  time  variable  t,  and  with  the  time-dependent  potential  V|  (x.t)  given  by 


V,(X,T) 


2  2  2" 
n  (x.x)  P  x  p 

2  ^"2". 


(7) 


For  the  ionosphere,  the  dielectric  constant  Is  given  by 


n2  (x.z)  =  e(x,z) 


vp(x-z) 

2' 

v 


(8) 


4 


where  vp  (x,z)  is  the  local  plasma  frequency  and  v  is  the  wave  frequency.  Methods  have  been 
developed  within  the  past  few  years  for  dealing  with  time-dependent  Hamiltonians  by  using 
the  method  of  Ermakov  Invariants2  or  equivalently  by  a  rescaling  of  space  and  time  variables 
in  conjunction  with  a  unitary  transformation  of  the  wave  function.2  4  The  most  general 
potential  that  can  be  dealt  with  by  these  methods  must  be  expressible  In  the  form 


Ve(X.T)  =  v  [(x  -  0(X)>  /P(T)]  /p2(t)  +  gQ(X)  +  X  .  g  ^X)  +  X2  •  g2(X). 


(9) 


Where  the  scaling  function  p(x)  and  shift  function  o(x)  must  obey  the  differential  constraint 
conditions 


P(X)  +  2g,(X)p(X)  a  K/p  (T) 


(10) 


o>)  +  2g2(x)o(x)+  g,(x)-0 


(11) 


with  tc  an  arbitrary  constant  and  with  go,  gi ,  and  g2  arbitrary  functions  of  time. 

The  equations  of  motion  for  a  ray  (or  a  classical  particle)  equivalent  to  the  wave  Eq.  (6) 
with  a  potential  given  by  Eq.  (9)  are 


d2x 

dt2 


'-2g2(T)x-gl(T)-^^>  -1 


3 

P  (*) 


(12) 


The  Ermakov  Invariant  for  this  system  Is  given  by 


I 


2 

I[(x-o)  p-p(x-o)]  +v((x-o)/p). 


(13) 


4  Lewis,  H.R.,  and  Leach,  P.G.L.  (1982)  Exact  Invariants  for  Time- Dependent  Nonlinear 
Hamiltonian  Systems.  Nonlinear  Problems:  Present  and  Future.  Bishop,  A.R.,  Campbell.  D.K.. 
and  Nicolaenko,  B.  (eds)  p.  133,  North  Holland. 
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If  x  and  t  are  transformed  as 


,  X  —  <3  (X) 
X  =  - — 

p(*) 


(14) 


and 


x 


dx" 

pV) 


(15) 


the  equations  of  motion  and  Ermakov  Invariant  are  transformed  to 


d2x'  ,  dv(x') 

-  s -KX  — - 

dt'2  dx' 


(16) 


and 


( dx' 
dr7 


+  v  (x') 


(17) 


respectively. 

We  will  call  the  (x'.x'l  coordinate  system  the  local  frame  and  (x,x)  the  global  frame.  From 
the  equation  of  motion  for  x'  given  by  Eq.  (16),  it  can  be  seen  that  in  the  local  frame  (for  the 
case  u  =  0)  the  scalar,  linear,  and  quadratic  portions  of  the  potential  proportional  to  gr  gj.  and 
g3  drop  out  and  that  only  the  potential  v  (x'),  which  can  be  thought  of  as  being  centered  about 
the  time  varying  position  of  its  minimum  x,^  *  o  (x),  affects  the  particle  motion.  The  local 
coordinate  x'  is  thus  simply  related  to  the  original  global  coordinate  x  by  a  time  varying  shift 
o  (t)  (which  locates  a  potential  extremum  or  minimum)  and  a  time  varying  scale  factor  p  (x). 

The  expression  for  the  Ermakov  invariant  in  terms  of  local  variables  given  by  Eq.  (17) 
appears  in  the  form  of  a  total  (local)  energy  that  is  conserved.  It  should  not  be  confused  with 
the  total  energy  of  the  system  given  by 
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E<,)M2(d7)2*V'(’‘’',) 


(18) 


which  Is  time  varying. 

For  the  wave  equation  the  Ermakov  Invariant  becomes  an  operator  with  the  replacements 


dx  _>  p 
dx 


(19) 


and 


dx'  ,  .  d 

->  p  =  - 1  — 

dx'  dx' 


(20) 


The  respective  Invariants 


I'  --1  —  +  1kx'2  +v(x') 
0,5  2  ax'2  2 


(21) 


and 


1 

op 


(22) 


are  related  by  a  unitary  transformation2 


I 

op 


e‘*I'  e-‘* 
op 


(23) 


where 


7 


<X>  =  &  +  xj  -  ox . 


(24) 


By  using  the  constraint  Eqs.  (10)  and  (11).  It  can  be  shown  that  both  Iop  and  I  are  truly 
Invariant,  that  Is. 


-d-1  =  I[l  . 

dt  °p  i  op 


Hi  +  i-i 
j  9t 


op 


(25) 


and 


Al  =  0. 

dt 


(26) 


3.  A  MORSE  POTENTIAL  MODEL  FOR  THE  E-F  LAYER  TROUGH 

In  order  to  apply  Ermakov  Invariant  theory  to  the  Ionospheric  ray  tracing  problem,  the 
ionospheric  potential  function  Vi  (x,t)  given  by  Eq.  (7)  must  be  cast  Into  the  form  of  the  lime 
dependent  function  VE  (x,i)  given  by  Eq.  (9).  Neglecting  for  the  moment  the  presence  of  an  Fj 
layer,  an  effective  ionospheric  potential  will  resemble  the  double-well  potential  VCH  (x.x) 
shown  In  Figure  1  that  was  formed  by  superposing  E  and  F2  layers  approximated  by  the 
Chapman  function 


fCH<X>  =  e 


l/2(l  -x-e'*) 


(27) 
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AVCH  M 


Figure  1.  Plot  of  the  Effective  Ionospheric  Potential 
VCH  (x)  as  a  Function  of  Height  x. 


In  terms  of  the  layer  half-widths  y,  (z),  yf  (z)  heights  xe  (z).  Xf  (z)  and  frequencies  \>e  (z).  (z) 

(provided  by  IONCAP)3  this  effective  Ionospheric  potential  Including  the  curvature  term  Is 
written  as 


(28) 


Finding  a  potential  function  approximating  Vch  (x.t)  with  known  analytic  solutions  and 
separately  adjustable  well  parameters  would  seem  to  be  a  remote  possibility.  If  the 
requirement  of  separately  adjustable  well  parameters  is  relaxed,  it  is,  however,  quite  possible 
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that  very  good  approximations  to  the  potential  of  Figure  1  could  be  found  by  employing  the 
Darboux  transformation.5  The  effect  of  the  Darboux  transformation  given  by 


V(x)  =  4>0(*) 


dx2 


(29) 


dx  dx  <j>o(X) 


(30) 


is  to  transform  a  single-well  potential  V  (x)  with  known  solution  $  (x)  Into  a  generally  non- 
symmetric  double-well  potential  V  with  solution  y  .  By  its  construction,  however,  the 
Darboux  transform  does  not  allow  either  well  to  vary  independently  of  the  other,  thus  limiting 
its  usefulness  for  ionospheric  ray  tracing  where  the  ground-E  layer  and  E-F2  layer  potential 
wells  do  vary  independently. 

In  order  to  insure  Independent  variation  of  the  potential  wells  (brought  about  by  the 
Independent  variation  of  the  constitutive  ionospheric  layers),  it  is  thus  necessary  to  piece 
together  truncated  wells  or  layers  along  common  boundaries.  In  this  section  we  will  look  at 
an  example  of  approximating  the  ionospheric  potential  by  piecing  together  two  potential  wells 
truncated  at  xr.  the  position  of  the  E-layer  peak;  while  in  Section  5  we  discuss  approximating 
the  Ionospheric  potential  by  piecing  together  either  Inverted  Morse  potentials  or  parabolic 
layers  connnected  by  straight  line  segments. 

While  piecing  together  two  truncated  potentials  leads  to  difficulties  In  generating  solutions 
to  the  wave  (Eq.  (6)]  for  states  near  the  top  of  the  barrier  separating  the  two  wells  due  to 
tunneling  and  mixing,  these  difficulties  do  not  exist  for  solutions  to  the  equation  of  motion 
(Eq.  (12)|  describing  ray  trajectories. 

The  steps  for  ionospheric  ray  tracing  employing  the  results  of  Ermakov  invariant  theory 
for  potential  wells  formed  by  Joining  together  separate  soluble  sectors  are  the  following: 

1)  Choose  a  form  VE  (x.t)  (Eq.  (9)1  that  can  closely  approximate  the  ionospheric  potential 
V,  (x.t)  (Eq.  (7)J  over  its  entire  range  by  varying  the  scale,  shift  and  level  functions  p  (t),  a  (t). 
and  gQ  (t)  appropriate  to  each  sector  with  g,  (t)  ,  and  gj  (t),  set  to  zero. 

2)  On  the  scale  of  a  wavelength  of  the  propagating  wave,  a  and  p  are  assumed  to  be 
negligible.  For  a  given  great  circle  path  p  (x|,  o  (x),  and  go  (x)  appropriate  to  each  sector  can  thus 

be  approximated  by  fitting  VE  (x,t)  to  V]  (x.t)  for  each  range  point  t,. 

3)  Launching  a  ray  at  some  point  (x,t0)  in  VE  (x,t)  at  a  given  incidence  angle  and  frequency 
determines  the  Ermakov  invariant  I  for  that  potential  sector  and,  assuming  the  solution  to  the 


5  Zheng.  W.M.  (1984)  The  Darboux  Tranformatlon  and  Solvable  Double- Well  Potential 
Models  for  Schroedlnger  Equations,  J.  Math.  Phys.,  25:68. 
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equation  of  motion  (Eq.  (16)1  (with  k  =*  0)  Is  known,  the  ray  path  Is  determined  up  to  the  point 
where  the  ray  crosses  a  boundary  into  another  sector  of  the  potential. 

4)  At  a  boundary  crossing  using  the  known  x.  x,  o,  p,  d,  and  p.  the  Ermakov  Invariant  Is 
calculated  for  the  new  sector  of  the  potential  and  motion  In  this  sector  Is  determined  as  In 
step  3.  Any  rays  that  hit  the  ground  at  x  =  0  are  assumed  to  be  perfectly  reflected. 

5)  The  final  ray  path  x  (x) .  Is  found  by  Inverting  Eqs.  (14)  and  (15)  and  piecing  the  ray 
paths  In  the  various  sectors  together. 

Instead  of  treating  the  entire  ionospheric  potential  at  once,  we  shall  first  focus  on  the 
potential  well  formed  by  the  E  and  F2  layers  and  the  propagation  of  rays  In  this  E-F2  duct. 
From  the  graph  of  VCH  (x.zj  In  Figure  1  It  can  be  seen  that  the  E-F2  potential  Is  asymmetric, 
with  a  slope  that  is  steeper  above  the  minimum  than  below  It.  A  soluble  potential  that  can  be 
adjusted  to  fit  the  E-F2  is  the  Morse  potential  given  by 


VM  <*'> 


2 

r 


(31) 


and  plotted  In  Figure  2. 
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Figure  2.  Plot  of  the  Morse  Potential  VM  (x)  =  y1  (e°“  -  l)2 
as  a  Function  of  Height  x  with  T2  *  1  and  a*  1 . 


6  Nieto,  M.M..  Simmons,  L.M.  (1979)  Coherent  States  for  General  Potentials  III.  Non¬ 
confining  One- Dimensional  Examples,  Phys.  Rev.  D,  20:1342. 


(34) 


A  (E')  -  (y2  •  E'/^/Oy2  -  E') 


and 


»C<E') 


W  (rJ-E-)l 


1/2 


(351 


For  E'  >  y2  the  solutions  are  given  by  the  same  expression  with  =  |  coc  |  replacing  toc  and 
cosh  replacing  cos. 

While  the  main  purpose  of  this  section  Is  to  model  and  examine  the  motion  of  rays  In  the 
E-F2  duct,  a  simplified  ground  E-layer  potential  will  be  added  on  In  order  to  have  a  rough 
model  available  for  launching  rays  from  a  ground  transmitter.  The  ground  E-layer  potential 


will  thus  be  approximated  by  a  simple  linear  potential  with  slope  g  ■  1/RQ  from  Vm  =  -  A  at 

2 

x  =  0  to  Vm  =  -  (1  +  ol/Ro)  at  the  minimum  x  «  ot  and  a  positive  slope  g  *  g,.  for  ol  i.  x  <  xe.  The 


right  side  of  this  potential  will  be  constrained  to  equal  the  value  of  VCh  (xe)  at  x  =  xe.  thus,  the 
position  of  the  minimum  oL  (*t)  Is  determined  by  the  Intersection  of  the  right  and  left  sides  of 
the  potential  and  will  move  lower  (higher)  as  the  value  of  VCh  (xe)  Increases  (decreases).  The 
general  solution  for  a  ray  In  this  potential  Is  given  by 


A-gt'2  +  v 
2  0 


x  +  x. 


(36) 


where 


{~  l/R|o  x'  <  0 


(37) 


To  cast  the  E-Fj  potential  well  Into  the  Ermakov  form,  the  Morse  potential  is  fit  at  each 
range  point  x(  by  a  least  squares  approximation.  It  Is  assumed  that  gt  (t),  g^  (x),  and  n  are  zero. 

At  the  Initial  range  point,  a  and  ^are  determined  by  least  squares  fit  with  p  (0)  =*  1.  go  (0)  »  0 
and  o  (0)  -  x,„|n.  At  subsequent  range  points,  o  (xj).  p  (x,).  and  gQ  (X|)  are  fit  by  minimizing 
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(*.',) -r2  (e-i-’WlW- if/ g0  (,,) 


dx 


(38) 


where  X|  and  xr  bound  the  E-F2  potential  as  described  Reference  7. 

Having  determined  the  parameters  a.  y3,  and  the  set  of  values  for  (o  (tt)|  and  (p  (xt)),  the 
path  for  a  ray  launched  in  the  E-F2  duct  at  some  arbitrary  point  (at  x  ■  xe  if  it  arrived  by 
passing  over  the  E-layer  peak)  will  be  given  for  all  subsequent  range  values  by  Eq.  (32)  until  it 
leaves  the  duct,  either  by  passing  over  the  E-layer  peak  at  x  »  xe  or  at  an  analagous  point  x  =  Xf 
which  is  taken  as  representing  the  F-layer  peak.  The  Ermakov  invariant  Im  determined  by  the 
values  of  the  quantities  x.  p.  and  d.  the  potential  VM  [(x  -  o)  /  pi  and  the  first  derivatives 
x.  p,  and  6  at  the  launch  point,  can  then  be  used  to  determine  the  velocity  at  the  exit  point  as 


x 


dx' 

dt' 


/  p  +  a  +  p  x' 


(39) 


where  dx'/dt'  is  related  to  IM  by 


1 


M 


+  VM(X'). 


(40) 


If  the  ray  crosses  over  into  the  ground-E  potential  well,  the  new  invariant 


.2 

,GS'y  +  [X  (%)]•«,■ 


(41) 


with  Tb  the  range  at  the  boundary  crossing,  is  conserved  as  long  as  the  ray  propagates  in  the 
ground  E-layer  potential.  If  we  were  able  to  treat  the  ionospheric  double- well  potential  exactly 
in  the  Ermakov  formalism,  there  would  be  only  one  invariant.  As  we  have  divided  the 
potential  into  sectors,  there  will  be  an  Ermakov  invariant  connected  with  propagation  in  each 
sector  and  crossing  from  one  to  the  other  will  entail  a  new  calculation  of  the  sector  Just 
entered.  Since  p  and  a  for  the  Morse  potential  can  change  Independently  of  p  (-  1)  and  a  for  the 
ground  E-layer  potential,  there  is  no  necessity  for  Im  (x  1 )  to  equal  Im  (t2)  if  the  ray  has 
propagated  in  the  ground  E-layer  potential  at  some  range  between  and  tj.  A  similar 
statement  holds  true  for  Iq  . 
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The  details  of  numerical  fitting  the  Morse  potential  wells  and  plotting  ray  paths  are 
discussed  in  Reference  7.  We  note  that  since  both  the  E  and  F2  layer  peaks  have  been 
approximated  by  functions  that  have  non-parabolic  slopes  at  xe  and  xf,  there  will  be  no 
Pedersen  high  rays. 


4.  COHERENT  8TATE8  AND  RAT  PATHS  POR  THE  MOR8E  POTENTIAL 


While  the  discussion  in  Section  3  was  concerned  with  determining  solutions  of  the 
classical  equations  of  motion  for  ray  paths,  the  same  time-dependent  theory  can  be  applied  to 
solutions  of  the  wave  [Eq.  (6)|  if  the  operator  form  for  the  Ermakov  invariants.  Eq.  (22)  or  (23) 
is  used.  Dhara  and  Lawande8  have  shown  that  the  Feyman  propagator  for  a  time-dependent 
system  that  obeys  the  constraint  conditions  (Eq.  (10)  and  (1 1)1  may  be  written  as 


MV  V  V ’.)  ■  Ov.)"1'2 «,k^-*»K0(,;.vb:  X-..,;) 


(42) 


where  KQ  is  the  propagator  in  the  local  potential  v  (x') 


Ko(*'b-'b:  *'.•«;) -x 


d-( 

• pmJ  -I  pw  J  ”1  p<\)  J 


(43) 


and  where  E'  are  the  eigenvalues  of  I '  with 

n  op 


X  ■ 


p*2 

2  P 


+  x  (d  -  op/p) 


if  P2('t')-i7(o/p)2dx'. 

2  Jo  dT 


(44) 


Since 


7  Yukon.  S.P.  (1986)  Calculation  of  Ray  Paths  in  the  Ionosphere  Using  an  Analytic  Ray 
Tracing  Technique.  RADC-TR-86-125,  ADA178888. 

8  Dhara.  A.K.,  and  Lawande.  S.V.  (1984)  Feynman  Propagators  for  Time  Dependent 
Lagrangians  Possessing  an  Invariant  Quadratic  In  Momentum.  J.  Phys.  A,  17:2423. 
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V  (v  %)  =  J  d  X.  K0  (V  V  X.  \)  *  (V  •'a)- 


(45) 


The  propagator  of  Eq.  (43)  projects  the  wave  function  y  (x..  tj  Into  a  sum  of  'local' 

eigenfunctions  6  (  —  each  of  which  has  Its  phase  Increased  according  to  Its  'local'  energy 
n\  p  ) 

eigenvalue  times  the  local  time 


t  = 


dt" 

P2(*f 


(46) 


If  an  appropriate  soluble  double-well  potential  could  be  found  that  would  model  the 
Ionospheric  potential  (Fj  layer  neglected),  then  the  above  propagator  would  convey  all  of  the 
Information  necessary  for  mapping  and  for  wave  propagation.  Since  It  Is  extremely  difficult 
to  find  a  soluble  double-well  potential  that  can  model  a  time-path-dependent  Ionosphere,  a 
practical  measure  Is  to  piece  together  truncated  soluble  single-well  potentials  as  outlined  In 
Section  3  The  resulting  total  wave  function  can  be  approximated  by  appropriate  super¬ 
positions  of  wave  functions  belonging  to  each  well.  If  eigenfunctions  and  eigenvalues  can  be 
found  for  each  truncated  well  separately,  the  combined  potential  wells  may  be  described  by  the 
same  formalism  that  is  used  for  treating  coupled  optical  waveguides.9 

If  we  can  confine  our  attention  to  low  lying  states  in  each  of  the  potential  wells  separately, 
the  details  of  how  the  potentials  are  truncated  should  not  greatly  affect  these  states.  For  low 
lying  states  in  the  E-Fa  potential  that  can  be  described  by  the  Morse  potential  (Eq.  (31)1  the 
eigenvalues  and  eigenfunctions  are  given  by10 


E 

n 


2k 


(n + 2)  -  (n + 2) 


(47) 


with 


9  Marcuse,  D.  (1974)  Theory  of  Dielectric  Optical  Waveguides.  Academic  Press,  New  York. 

10  Morse.  P.M.  (1929)  Diatomic  Molecules  According  to  the  Wave  Mechanics  D. 
Vibrational  Levels.  Phys.  Rev.,  34:57. 
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(48) 


E0  =  h2  a2/ 2  m ,  X2  =  Y2/EQ 


and 


*n  -  N  (n,  X)  4°  <y) 


(49) 


a  (2  X  2  n  1 )  p  ✓  + 
r(2X-n) 


1/2 


y(X-  1/2-n)  e-y/2 


^(2  V-2n 
n 


i) 


<y) 


with 


y  =  2  X  e~  “  x 


(50) 


where  L(m) 
n 


Is  a  generalized  Laguerre  polynomial. 


4.1  Coherent  8tates  of  the  Morse  Potential 

It  has  been  pointed  out  by  Schrodlnger11  and  subsequently  by  several  other  authors12  that 
for  wave  packets  formed  to  minimize  the  position-momentum  uncertainty  for  a  harmonic 
oscillator,  the  expectation  value  <x  (x)>  will  follow  the  classical  motion  of  an  harmonic 
oscillator.  Subsequent  authors  have  generalized  such  minimum  uncertainty  states  to  coherent 
states  that  can  be  defined  either  as  eigenstates  of  the  annihilation  operator 


a 


(p  +  IBflt  xy  V2m© 


(51) 


11  Schrodlnger,  E.  (1926)  NaturwissenschaAen.  14:664.  English  translation  In  Schrodlnger. 
E.  (1928)  Collected  Papers  on  Wave  Mechanics,  p.  41,  Blackle  &  Son,  London. 

12  Nieto,  M.M.,  Simmons  Jr.,  L.M.  (1979)  Coherent  States  for  General  Potential?,  I. 
Formalism,  Phys.  Rev.,  D20:1321. 
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such  that  a  |  a  >  =  a  |  a  >.  or  as  a  displacement  by  a  of  the  ground  state  1 0  >. 

For  the  harmonic  oscillator  the  two  definitions  are  equivalent  due  to  the  special 
circumstance  that  the  oscillator  frequency  Is  Independent  of  energy,  with  elgenenergles 
equally  spaced.  For  the  Morse  potential  these  conditions  do  not  hold  as  can  be  seen  from  the 
expression  for  the  frequency  [Eq.  (35)1  which  is  energy  dependent  and  from  the  expression 
[Eq.  (47)]  for  En  which  depends  on  n2. 

By  treating 


X  =  c~ax  -  --Y  -  -  A  (E)  sin  [«c  (E)  t] 

ly  -  e) 


(52) 


and 


Pc  =  aPe  aK  =  m  coc  (E)  cos  [a^  (E)  t] 


(53) 


as  an  appropriate  pair  of  classical  (but  not  canonical)  variables  for  the  problem,  since  their 
classical  equations  of  motion  are  the  same  as  the  harmonic  oscillator.  Nieto12  formed 
minimum- uncertainty  coherent  states  that  minimize  <  A  Xc  ><  A  Pc  >.  In  tracing  out  this 
subsequent  time  evolution  of  these  states  it  was  found  that  after  a  few  cycles  the  wave  packets 
would  tend  to  break  up  and  disperse.  More  recently  Gerry13  has  shown  that  the  Morse 
oscillator  may  be  treated  as  an  harmonic  oscillator  evolving  In  a  curved  phase  space  with  a 
transformed  time  parameter 


x 


dt'e" 


■(*') 


(54) 


The  coherent  states  can  thus  be  formed  as  minimum  uncertainty  states  In  the  transformed 
harmonic  oscillator  representation.  When  projected  back  into  original  variables,  <  x  (x)  >  for 
these  coherent  states  should  follow  the  classical  path  for  a  Morse  oscillator. 


13  Gerry,  C.C.  (1986)  Coherent  States  and  a  Path  Integral  for  the  Morse  Oscillator,  Phys. 
Rev.,  A33:2207. 
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4.2  Focusing  Properties  of  the  Morse  Potential 


It  has  been  noted  for  many  years14  that  underwater  sound  channels  exist  In  the  oceans 
that  can  play  the  same  role  for  acoustic  waves  as  Ionospheric  ducts  do  for  HF  electromagnetic 
waves.18  In  the  underwater  sound  case  there  exists  a  focusing  phenomenon  shown  In  Figure  3 
wherein  rays  launched  from  a  source  near  the  sound  axis  will  be  repeatedly  focused  as  they 
travel  downrange  In  the  sound  channel.  The  regions  where  this  quasi-focusing  occurs  are 
called  convergence  zones.  The  sound  speed  profile  for  underwater  sound  channels  Is  well- 
approximated  by  the  Munk  profile16 


X  -  x_ 


(55) 


shown  In  Figure  3.  Near  the  origin  the  Munk  profile  can  be  seen  to  be  quite  similar  to  the 
Morse  potential  for  small  q.  Writing  the  Morse  potential  as 


VM  -  y2  (e11  -  2el'/1+  1 ) 


(56) 


and  expanding  only  the  second  term  to  first  order,  that  Is 


-2en/2  =  -2-q 


(57) 


yields  an  expression  which  Is  Identical  to  the  Munk  profile.  Repeated  focusing  could  therefore 
be  expected  to  occur  In  the  case  of  ducted  Ionospheric  propagation  described  by  the  Morse 
potential. 


14  Ewing.  W.M.,  Worzel,  J.L.  (1948)  Long  Range  Sound  Transmission,  GeoL  Soc.  Am 
Mem.  27.  Part  m.  1. 

18  Quack.  E.  (1927)  Propagation  of  Short  Waves  Around  the  Earth,  Proc.  IRE,  15:341.  Further 
references  In  Toman,  K.  (1979)  High-Frequency  Ionospheric  Ducting  -  A  Review,  Radio  Science. 
14:447. 

16  Brekhou sklkh ,  L.  and  Lyzanov,  Yu  (1982)  Fundamentals  of  Ocean  Acoustics,  Springer- 
Verlag.  Berlin,  Heidelberg  New  York,  p.  118. 
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Co  Cft 


Zones  from  Reference  16. 


As  It  turns  out.  repeated  focusing  In  the  Morse  potential  can  be  exact,  that  Is  real  point 
focusing  as  opposed  to  progressively  wider  convergence  zones.  To  see  this  we  first  note  that  to 
have  repeated  focusing,  the  propagation  constants  pn  must  satisfy 


P  -  P  ,  =  I 

rn  rn  nn  l 


(58) 


where  I  ,  Is  an  arbitrary  Integer  function  and  L  Is  the  focal  length.17  There  exist  families  of 
refraction  Index  profiles  that  are  finite  at  spatial  Infinity  and  have  the  focusing  property 
described  by  Eq.  (58),  for  example 


2 

Vf (x)  =  sech  (A(x-xq))  where  pn  =  n»k0. 


(59) 


For  the  Morse  potential,  which  does  not  fit  this  category  since  VM  (x)  Is  unbounded  as 
x  ->  +  It  Is  still  possible  to  satisfy  the  focusing  condition  [Eq.  (58)]  due  to  the  presence  of  the 


17  Yukon.  S.P.  and  Bendow.  B.  (1980)  Design  of  Waveguides  with  Prescribed  Propagation 
Constants.  J.  Opt  Soc.  Am.,  70:172. 
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(n  +  1/2 )2  term  In  the  expression  for  the  energy  eigenvalue  given  In  Eq.  (47).  In  order  to  satsify 
condition  1^-  =  1 .  (3n  can  be  written  as 

Pn  =  p  +  q»n.  (60) 


This  relation  will  then  hold  if  /  =•  /2.  where  n*/2  is  the  magnitude  of  the  potential  well 

minimum.  From  the  ionospheric  paths  that  so  far  have  been  fitted  by  the  Morse  potential  it 
appears  that  this  condition  can  be  met,  but  only  for  frequencies  in  the  range  of  6-10  MHz 

Since  only  a  portion  of  the  E-F2  Ionospheric  duct  can  be  modeled  by  the  Morse  potential, 
repeated  focusing  cannot  be  expected  to  hold  for  weakly  bound  ducted  modes,  but  should  apply 
to  the  tightly  bound  modes.  From  the  correspondence  between  rays  and  coherent  states 
discussed  In  the  previous  section  It  is  expected  that  repeated  focusing  will  thus  hold  for  low 
angle  ducted  rays. 

It  can  also  be  seen  that  repeated  focusing  will  be  an  invariant  j  ,-ertv  of  an  ionospheric 
duct;  If  repeated  focusing  holds  true  In  the  local  frame,  a  som*  vhat  distorted  repeated  focusing 
will  hold  in  the  global  frame.  Thus,  rays  leaving  n  initial  point 


xo”°(to)+  xo*P(To) 


(61) 


at  range  to  will  refocus  at 


x,  =°(Ti)+  x'i'P(T1) 


(62) 


at  range  ii,  where  t0  and  tj  are  related  by 


L  * 


dt 

P2(t) 


(63) 


providing  the  differential  constraint  conditions  (Eq.  (10)  and  Eq.  (11)|  are  satisfied.  In  the 
global  frame  the  focal  length  Is  a  function  of  p  (r),  and  the  mapping  of  an  Interval  a  <  x  s  b  will 
be  displaced  by  a  (x  i)  -  a  (t0)  and  enlarged  or  reduced  depending  on  whether  p  (tj)  Is  greater  than 
or  less  than  p  (x2)  . 
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5.  RAY  TRACING  IN  THE  IONCAP  IONOSPHERE 


The  Ionospheric  electron  density  profiles  provided  by  the  IONCAP  program3  consist  of  E. 
F2  and  an  optional  Fi  layers,  Joined  together  by  straight  line  segments.  The  parabolic  layers 
are  described  by  the  plasma  frequency 


up(x)='>o[1-(x-xo)2/y?] 


(64) 


where  \>o.  xo.  and  yo  are  the  peak  plasma  frequency,  layer  height  and  half  width  respectively 

2 

with  the  refractive  Index  e  (x)  =  n2  (x)  related  to  \>p  (x)  by 


«  (*)  -  1  -  \  (x)  A>2- 


(65) 


The  base  of  the  E-layer  Is  set  at  hj  =  90  km.  The  straight  line  that  models  the  valley  between 
the  E-  and  F2-  layers,  starts  at  a  point  h2  where  the  plasma  frequency  Is  0.8516  £de  and  ends  at 
a  point  h3  where  the  plasma  frequency  Is  0.98  foE-  'This  yields 


=  120.48  km 


(66) 


and 


h 


3 


xof2~  y0F2 


VT 


(67) 


The  effective  ionospheric  potential  is  that  given  by  the  sum  of  curvature  correction  term 


5V  -_L 
E  2 


( 


+  2x_lx_ 

2  R2 
0  Roy 


(68) 
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and 


=  -^  (°2  (x.t)  +  p2) 


(69) 


where  p  =  cos  (0).  with  0  the  launch  angle  with  respect  to  the  horizontal  direction.  The 
resulting  potential  well  Is  shown  In  Figure  4.  If  the  Ionospheric  potential  is  thought  of  as  two 
separate  wells  Joined  at  the  E-layer  peak  as  In  the  case  of  the  Morse  oscillator  model  of 
Section  3.  each  separate  well  Is  then  composed  of  parabolic  and  linear  pieces,  and  would  not  be 
in  the  class  of  potentials  soluble  by  Ermakov  theory  since  the  shape  of  the  well  cannot  be 
expressed  as  a  single  function  V  |(x  -  a)  /  p]. 


Figure  4.  Plot  of  the  IONCAP  Effective  Ionospheric 
Potential  VE  (x)  as  a  Function  of  Height  x. 
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If.  however,  the  Ionospheric  potential  Is  broken  up  into  sectors  as  shown  tn  Figure  5. 
Ermakov  theory  may  then  be  applied  to  each  of  the  four  potentials.  The  potentials  in  Sectors 
II  and  IV  are  not  what  one  usually  thinks  of  as  potentials  since  they  have  no  local  minimum, 
but  they  are  valid  potentials  for  Ermakov  theory  and  have  recently  been  shown18  to  have 
coherent  states  that  follow  the  classical  motion.  For  an  Inverted  parabola  centered  about  x  =  0 
the  classical  solutions  are 


x  (x)  =  x0cos  (tot)  +  pQ  sin  (tot) 


(70) 


for  E=  1/2  (pj-xj)<0 


and 


x  (t)  =  xQ  ch  (to  t)  +  pQ  sh  (tot) 


(71) 


for  E  >  0. 


10  Barton.  G.  (1986)  Quantum  Mechanics  of  the  Inverted  Oscillator  Potential.  Ann.  ofPhys., 
166:322. 
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Figure  5.  The  lONCAP  Potential  of  Figure  4  Broken 
up  Into  Sector  Potential  Wells  I  through  IV. 


Since  \>o.  x0,  and  yo  are  provided  by  IONCAP  at  each  range  point,  o  (z).  and  p  (z)  can  be 
found  algebraically  and  need  not  be  determined  by  the  fitting  procedure  that  was  necessary  for 
the  Morse  potential.  The  tradeoff  for  not  having  to  fit  the  parameters  a  (z).  and  p  (z)  Is  that 
there  are  now  more  boundary  crossings  (three  vs  one  for  the  Morse  case)  where  the  portion  of 
the  range  step  that  Is  In  each  sector  must  be  determined;  for  example  for  x  (zt)  In  region  III  and 
x  (zi  +  1)  In  region  IV  It  Is  necessary  to  find  zt,  such  that  x  (z^  =  h3  (zb).  Determining  Zb  Is 
straightforward  for  the  above  example  (complicated  somewhat  by  the  fact  that  the  h3  boundary 
is  a  function  of  range).  If  the  ray  did  not  reverse  direction  between  zt  and  z(  +  1.  If  there  Is  a 
velocity  reversal  In  this  interval  then  three  possible  cases  must  be  considered  as  depicted  In 
Figure  6  below. 
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Ill 


IV 


Figure  6.  Three  Possible  Cases  for  Turning  Point 
Behavior  Near  a  Sector  Boundary. 


To  distinguish  between  them  it  is  thus  necessary  to  calculate  the  turning  point. 

One  further  situation  where  care  must  be  exercised  is  for  rays  that  cross  very  near  the  peak 
of  the  nigh  time  E-layer.  The  peak  of  the  parabola  representing  the  E- layer  electron  density 
will,  upon  being  added  to  the  curvature  term  to  form  the  effective  potential,  be  shifted  to  a 
lower  altitude  by  an  amount 


as  well  as  being  shifted  down  In  magnitude  by  an  amount 
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(73) 


It  is  possible  for  the  potential  peak  to  be  shifted  down  (In  altitude)  so  far  that  the  effective 
potential  between  h  i  and  I12  is  the  mapping  of  a  topside  segment  of  the  inverted  parabolic 
potential  and  does  not  contain  the  peak.  A  ray  that  crosses  from  x  <  hj  to  x  >  hi  may  then 
find  Itself  on  the  topside  region  of  an  Inverted  parabolic  potential  even  though  its  energy  is 
negative  with  respect  to  the  shifted  peak.  While  unexpected,  such  behavior  is  easily  dealt  with 
simply  by  checking  whether 


h;  =  (hi-o(Tb))/p(Tb) 


(74) 


is  greater  than  or  less  than  zero.  For  h  j  <  0  the  E  layer  boundary  is  still  left  of  the  peak  in 
local  coordinates,  while  for  hj  >0  the  ray  is  on  the  topside  of  an  inverted  parabolic  potential. 

The  testing  to  determine  initial  conditions  occurs  only  for  range  steps  where  the  ray 
crosses  a  boundary.  Thus  given  the  ionospheric  parameters  for  a  given  great  circle  path  from 
IONCAP.  tracing  rays  employing  Ermakov  invariants  essentially  amounts  to  Joining  up  the 
analytic  solutions  In  each  sector.  The  only  real  numerical  work  involved  is  in  the  calculation 
of  the  fractional  range  step  remainders  in  crossing  sector  boundaries.  For  an  around-the- 
world  path  there  can  be  on  the  order  of  20  complete  oscillations  In  the  EF  potential  well  which 
would  require  40  numerical  minimizations  to  determine  the  range  step  remainders.  These 
may  be  accomplished  very  efficiently  using  standard  Newton’s  minimization  methods  so  that 
a  complete  round-the-world  ray  trace  may  be  carried  out  in  approximately  one-half  second. 

An  tonogram  for  the  same  path  covering  angles  from  0°  to  40°  In  one-half  degree  steps  and  the 
even  frequencies  from  6  to  30  MHz  takes  approximately  3  minutes  to  calculate. 

The  F  i  layer  as  modeled  by  IONCAP  may  also  be  straightforwardly  Included  in  the  ray 
tracing  code  as  a  third  parabolic  (or  linear)  layer  centered  at  the  F]  layer  height.  This  will 
entail  the  calculation  of  one  more  boundary  crossing  for  motion  between  the  E  and  F2  layers. 
Since  a  ray  may  become  trapped  between  the  F  \  and  F2  boundary,  the  same  checking  for  one- 
step  ray  reversals  at  the  F1-F2  boundary  must  be  carried  out  as  was  done  for  the  case  of  the 
III-F2  boundary  described  above.  Thus  an  approximate  10  percent-30  percent  increase  in 
computation  time  could  be  expected  by  turning  on  the  Fj  layer. 

Some  examples  of  ray  tracings  and  lonograms  made  using  the  methods  outlined  above  are 
shown  In  Figures  7  through  1 1  along  with  results  obtained  using  the  Jones-Stephenson 
program  for  the  same  IONCAP  parameters. 
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Figure  7.  Ray  Traces  Using  Current  Model  vs  Jones-Stephenson's1 
Program  for  Frequency  ■  16  MHz,  Takeoff  Angle  -  7.1°. 
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Figure  8.  Ray  Trace*  Using  Current  Model  vs  Jones- Stephen  son’s* 
Program  for  Frequency  -  16  MHz,  Takeoff  Angle  ■  6.7°. 
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Figure  9.  Ray  Traces  Using  Current  Model  vs  Jones-Stephenson's 1 
Program  for  Frequency  ■  16  MHz.  Takeoff  Angle  ■  6.4°. 
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SEASON  -  MARCH  EQUINOX  DIRECTION  -  DUE  NORTH 
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Figure  10.  lonogram  Produced  Using  Current  Ray  Tracing  Model 
with  Transmitter  Location  and  Time  as  Shown  In  Figure  7. 
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TRANSMITTER  -  ALICE  SPRINGS,  AUSTRALIA  TIME- 0300  UT 

SEASON  -  MARCH  EQUINOX  DIRECTION  -  DUE  NORTH 


Figure  11.  Ionogram  Produced  Using  Jones- Stephenson1  Program 
with  Transmitter  Location  and  Time  as  Shown  In  Figure  7. 
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Rome  Air  Development  Center 


RADC  plans  and  executes  research,  development,  test  and 
selected  acquisition  programs  in  support  of  Command,  Control, 
Communications  and  Intelligence  (CiI)  activities.  Technical  and 
engineering  support  within  areas  of  competence  is  provided  to 
ESD  Program  Offices  (POs)  and  other  ESD  elements  to 
perform  effective  acquisition  of  CSI  systems.  The  areas  of 
technical  competence  include  communications,  command  and 
control,  battle  management  information  processing,  surveillance 
sensors,  intelligence  data  collection  and  handling,  solid  state 
sciences,  electromagnetics,  and  propagation,  and  electronic 
reliability /maintainability  and  compatibility. 


